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ABSTRACT 

Seven years of pulse time-of-arrival measurements have been coUected from 
observations of the young pulsar PSR B2334+61 using the Nanshan radio tele- 
scope of Urumqi Observatory. A phase-connected solution has been obtained 
over the whole data span, 2002 August to 2009 August. This includes a very 
large glitch that occurred between 2005 August 26 and September 8 (MJDs 
53608 and 53621). The relative increase in rotational frequency for this glitch, 
I^Vgjv ~ 20.5 X 10~^, is the largest ever seen. Although accounting for less 
than 1% of the glitch, there were two well-defined exponential decay terms with 
time constants of 21 and 147 days respectively. There was also a large long-term 
increase in the spindown rate with Az>p/z> ~ 0.011 at the time of the glitch. A 
highly significant oscillation with a period of close to one year is seen in the post- 
glitch residuals. It is very unlikely that this can be accounted for by a pulsar 
position error or proper motion - it appears to result from effects interior to the 
neutron star. Implications of these results for pulsar glitch models are discussed. 

Subject headings: stars: neutron — pulsars: general — pulsars: individual (PSR 
B2334+61) 

1. Introduction 

Pulsars are well known for their very stable spin periods. However, long-term timing 
observations of young pulsars have revealed two categories of rotational irregularities: timing 
noise which is usually characterized by continuous and erratic fiuctuation in rotation period, 
and glitches which are sudden increases in their rotation rates. Typical increases in rotational 



frequency during a glitch are of the order of lS.Vglv ^ 10~^ to 10~^. Large glitches are 
predominately found in young pulsars with characteristic ages between 10 kyr and 1000 kyr 
and are often followed by a relaxation process in which some of the glitch decays. It is believed 
that glitches and post-glitch variations reflect the dynamics of the interior of the neutron star 
rather than magnetospheric phenomena. Observing glitches and measuring their subsequent 
decay processes, in addition to providing insights into the phenomenon itself, offer one of 
the few probes of the neutron star structure and thus the physics of ultra-dense matter (e.g., 
Larson & Link 2002; Alpar & Baykal 2006). 

The young pulsar PSR B2334+61 (PSR J2337+6151) was discovered in the Princeton- 
NRAO survey using the 92-m radio telescope at Green Bank (Dewey et al. 1985); it is 
located near the center of, and probably associated with, the supernova remnant Gl 14.3+0.3 
(Fiirst et al. 1993). The pulsar has a rotational period P ~ 0.495 s and a period derivative 
P ~ 1.9x 10"-*^^, implying a characteristic age r^ = P/{2P) ~ 41 kyr. B2334+6rs spin-down 
power, E = — 47r^/yz>, is 6.23x10^^ erg s"-*^, where the moment of inertia / is taken to be 
]^q45 g pj^2_ ^i^jg suggests that the pulsar may be a 7-ray emitter, although its distance of 
^ 3 kpc puts it at the bottom end of the E^^'^dr'^ metric for pulsars detectable by the Fermi 
Space Telescope (Abdo et al. 2010). X-ray emission from the pulsar was first detected in a 
ROSAT pointing (Becker et al. 1996) and was shown to be thermal using an XMM-Newton 
observation (McGowan et al. 2006). 

In this Letter we present the timing solution of PSR B2334+61 obtained from seven 
years of observation with the Nanshan 25-m telescope. The timing analysis reveals the largest 
known glitch ever observed, with a fractional frequency increase of APg/i' ~ 20.5 x 10~^, a 
double exponential decay of a small part of the glitch, a large permanent increase in spindown 
rate and a quasi-sinusoidal oscillation apparently triggered by the glitch. We give a short 
discussion about the implications of these results for glitch theories. 



2. OBSERVATION AND ANALYSIS 

Timing observations of PSR B2334+61 using the 25-m Nanshan radio telescope at 
Urumqi Observatory commenced in 2002 July. The observing system is described in Wang 
et al. (2001). In brief, the two hands of circular polarization at a frequency of 1540 MHz are 
fed through a 2 x 128 x 2.5 MHz-channel filter bank and digitized at 1 ms intervals. The 
observation time for this pulsar was 16 minutes. 



Off-line data reduction was performed in stages using the PSRCHIVE package^. The 
data were dedispersed and summed to produce a total intensity profile. Local times-of- 
arrival (TOAs) were determined by correlating this profile with a standard pulse profile of 
high signal-to-noise ratio. Correction of TOAs to the Solar-System barycentre and fitting 
of the timing model were done using TEMP02^ with the Jet Propulsion Laboratory DE405 
ephemeris (Standish 1998). TOAs were weighted by the inverse-square of their estimated 
uncertainty. The basic timing model is the Taylor series 

4>{t) = 00 + ^{t - to) + ]^p{t - hf + ^i){t - hf + ^•z>-(i - hf (1) 

where 0o is the pulsar phase at the reference bary centric time to. Glitches can be modelled 
as combinations of changes in z/, z> and i>, of which parts may recover on various timescales: 

v{t) = uo{t) + Aup + Aupt + -Aupf + Aud e'^l^^ (2) 

v{t) = Po{t) + APp + Ai)pt + Avd e-^l^^ (3) 

v{t) = uo{t) + Aup + Aud e-^l^^ (4) 

where Avp and Avp are permanent changes in v and v relative to the pre-glitch solution, Avd 
is the amplitude of a decaying component with a time constant of r^, Av = —AvdJTd and 
Avd = Afd/Td- The total frequency jump at the time of the glitch is Aug = APp + Afd- The 
degree of recovery is often described by the parameter Q = Aud/Aug. If required, multiple 
decay terms with different decay time constants can be fitted; the total frequency jump is 
then the sum of the permanent and decaying frequency jumps, and similarly for z> and u. 



3. RESULTS 

The data presented here were recorded between 2002 July 31 and 2009 August 1 (MJDs 
52486 and 55044). Figure 1 gives an overview of the results of a timing analysis, showing 
that a large glitch with Avg Ri 41 x 10~^ Hz occurred between 2005 August 26 and September 
7 (MJDs 53608 and 53621). The regular timing observations of this pulsar at Nanshan allow 
us to track the recovery process in detail. Most of the frequency jump Aug persists beyond 
the end of the data span, but the expanded Au plot given in Figure 1(b) and the plot of 
z> given in Figure 1(c) show that there was an initial exponential decay with a timescale of 



■^See http://psrchive.sourceforge.net 

^See http://www.atnf.csiro.au/research/pulsar/tempo2 



^ 150 days and an apparently permanent increase in slow-down rate \u\. As is described 
in detail below, there was an additional more rapid exponential decay with a timescale of 
^ 21 days following the glitch. 

The rotational and positional parameters for PSR B2334+61 from fitting of the timing 
model, Equation 1, to the pre- and post-glitch data are given in Table 1. Pulsar frequencies 
are in TDB units and uncertainties given in parentheses are in the last quoted digit and 
are twice the values given by tempo2. The position was derived from the pre-glitch data 
assuming the proper motion given by Hobbs et al. (2004) and held fixed for the post-glitch 
solution. To avoid most of the short-term post-glitch decay, the post-glitch solution was 
fitted to data starting approximately one year after the glitch. Even with this, there were 
significant systematic period variations beyond the second-derivative term, most of which 
were absorbed by a third-derivative term. 

Figure 1(c) shows that there was a very significant long-term increase in \u\ at the time 
of the glitch; the values of z> given in Table 1 refiect this change. Furthermore, there was 
a significant long-term change in i) at the time of the glitch, with the rate of decrease of 
\u\ increasing by more than four after the glitch. The corresponding braking indices are 
10.5 ± 0.2 and 46.8 ± 0.3 before and after the glitch, respectively. 

Table 2 gives glitch parameters from a fitting of a model to the timing data. Because 
of the large amplitude of the glitch, the final fit was obtained in several stages. For the 
first stage, fits the data span were terminated at MJD 54000, just over a year after the 
glitch, to avoid contamination by the long-term decay in z>. Az>p was held fixed at the 
value given by subtraction of the pre-glitch value of z> given in Table 1 from the post-glitch 
value, viz, —7.625 x lO"-*^^ s~^. Then, based on Figure 1, a single decay with an assumed 
time constant of 150 d was fitted to the TOAs, giving values of Aup ^ 4.1203 x 10~^ Hz 
and Audi ~ 1-30 x 10"'' Hz. Post-fit residuals from this fit are shown in the upper part 
of Figure 2. It is clear from this plot that there was an additional initial decay with time 
constant ^ 20 d. Fitting for this with the parameters of the longer decay held fixed gave 
Ah'd2 ^ 1-9 X 10~^ Hz and Td2 = 18.5 d. Residuals for this fit shown in the lower part of 
Figure 2 demonstrate that the observed TOAs over this data span are very well modelled 
by increments in u and z> together with two exponential decays with time constants around 
20 d and 150 d respectively. 

Using these parameters as starting points, it was then possible to fit the entire dataset 
from 2002 to 2009 to give the glitch parameters listed in Table 2. The pulsar position, 
proper motion, z^, z> and i> values were held fixed at the pre-glitch values from Table 1. 
All other parameters were solved for with the fit being iterated several times to ensure 



convergence^. There was a gap of 12 days between observations around the gHtch; 
the assumed gHtch epoch is in the middle of this gap. Although most derived 
parameters are essentially independent of the assumed epoch, the amplitudes of 
the decaying components do change significantly. Quoted uncertainties in the 
glitch parameters include a contribution from the uncertainty in the glitch epoch 
Residuals from this final fit are shown in Figure 3. Values for Az>p and Aup are close to the 
diflferences between z> and i) for the post- and pre-glitch fits (Table 1) as expected. They 
do not coincide exactly (within the estimated errors) mainly because the exponential terms 
included in the final fit absorb some of these terms. 

This final fit clearly represents the glitch and the subsequent decay extremely well. To 
put the scale of the residuals in Figure 3 into context, Aug represents a phase gradient of one 
turn every 6.2 hours and the Az>p and Ai>p terms correspond to a total phase contribution 
over the post-glitch data span of 66.2 and 2.76 turns respectively. The final rms residual 
corresponds to 0.0035 turns. Although good, the final fit is not perfect, with a reduced 
X^ of 3.32. The x^ is dominated by systematic deviations from the fit after the glitch. 
These appear to be quasi-sinusoidal; a Lomb-Scargle spectral analysis shows a periodicity 
in the post-glitch post-fit residuals of 364 d with an uncertainty of ~ 5 d and a significance 
(false-alarm probability) of 0.0005. The fact that the period of the post-glitch oscillations is 
equal to a year within the uncertainty immediately suggests that the oscillations are due to a 
position error or proper motion. However, we believe that this is not the correct explanation. 

The position derived from the pre-glitch data (Table 1) is fully consistent with the 
position extrapolated to MJD 53100 based on the position and proper motion obtained 
by Hobbs et al. (2004) using Jodrell Bank data spanning 16 years: 23*^ 37™ 05':779(28), 
+61°51'01'.'54(17). A fit of position only to the post-glitch data, holding all glitch and pulse 
frequency parameters constant at the final-fit values, moves the position by +O!059 of right 
ascension and +0'.''36 of declination, about four times the combined uncertainties. Compared 
to the extrapolated Hobbs et al. (2004) position, the shifts are about twice the combined 
uncertainties. Including position, proper motion, z^, z>, u and all glitch terms in a fit to the 
whole data set produced small changes in all glitch parameters (less than the uncertainties 
given in Table 2) and gave a proper motion of +111(24) mas yr"-*^ in the right ascension 
direction and +75(23) mas yr"-*^ in declination. These values are inconsistent with the Hobbs 
et al. (2004) values. Furthermore, the final post-fit residuals have significant annual terms 
of opposite sign before and after the glitch. These considerations strongly suggest that the 
solution illustrated in Figure 3 is the best possible, that the oscillations are confined to 



^Solving for the pulsar position, v, v and i) also converged well and did not change any of the parameters 
outside their uncertainty range, demonstrating that the covariance between the various parameters is small. 



post-glitch epochs and that they are not caused by a position error or proper motion of the 
pulsar. 

There is no evidence for a change in pulse shape at the time of the glitch, but the limits 
are rather poor as the typical signal-to-noise ratio from each 16-minute observation was only 
about 10. 



4. DISCUSSION 

The glitch in PSR B2334+61 is the largest in Az^g/z^ terms (and equal largest in /S.Pg 
terms) ever observed for any radio pulsar and probably in any pulsar (including magnetars).'^ 
With a characteristic age of 41 kyr it is in the age range of pulsars known to exhibit large 
glitches (see, e.g.. Yuan et al. 2010). In fact, it lies very close in the P — P plane to 
PSR J1806-2125 which in 1998 had a glitch with /\vg/p - 16 x 10"^ (Hobbs et al. 2002). 
There are other six pulsars close to PSR B2334+61 in the P — P plane: PSRs J1551— 5310, 
J1737-3137, B1758-23, J1838-0453, J1841-0524 and B2000+32. For PSR B1758-23, six 
glitches of moderate size (Az^^/z^ in the range (0.017 — 0.35) x 10~^) have been detected. 
For PSRs J1737-3137 and J1841-0524 a single glitch of relative size - 1.3 x 10"^ and 
1.0 X 10~^ respectively have been reported (Weltevrede et al. 2010). Other pulsars have not 
been reported to glitch so far. This is the first glitch reported for PSR B2334+61; since the 
pulsar was observed for 16 years from 1987 by Hobbs et al. (2004), the time between glitches 
is at least 22 years. 

In contrast to the case for PSR J1806— 2125, for PSR B2334+61 observations were 
obtained close to the date of the glitch allowing a detailed examination of the post-glitch 
behaviour. PSR B2334+61 continues the trend for minimal post-glitch recovery of Az^g, 
i.e., small Q, for large glitches in middle-aged pulsars (Lyne et al. 2000; Yuan et al. 2010). 
Despite this, the good sampling of the post-glitch behaviour shows clear evidence for two 
distinct exponential recoveries and, together with permanent jumps in z> and i>, these provide 
a remarkably accurate description of the post-glitch behaviour. Previously, only the Vela 
pulsar has shown clear evidence for multiple exponential decays associated with a given glitch 
(Dodson et al. 2002). For Vela, the longest decay timescale observed by Dodson et al. (2002) 
is about 19 d, comparable to the shorter decay timescale observed for PSR B2334+61. There 
is evidence for longer exponential decays in the Vela pulsar glitches which ultimately become 
more linear (Lyne et al. 1996; Lyne et al. 2000; Wang et al. 2000). 



*See the Glitch Table in the ATNF Pulsar Catalogue (www.atnf.csiro.au/research/pulsar/psrcat). 



The standard model for large pulsar glitches envisages a sudden transfer of angular 
momentum from a superfluid component to the rest of the star - the vortex unpinning 
model (Anderson & Itoh 1975; Alpar et al. 1984). The theoretical understanding of glitches 
includes the trigger mechanism for the glitch, the strength of the vortex pinning and the 
post-glitch evolution. Larson & Link (2002) investigated a model in which vortex unpinning 
results from a sudden heating of the star, possibly due to a starquake or crustal movement, 
whereas Glampedakis & Andersson (2009) suggested that an r-mode instability may trigger 
a global unpinning of vortices leading to a glitch. Alpar et al. (1993) and earlier papers by 
these authors attribute the long-term post-glitch relaxation to vortex creep, with different 
regions of the star having different properties to account for the sometimes complicated post- 
glitch behaviour. Ruderman et al. (1998) suggested that crustal motions could account for 
essentially all glitch properties; see also Ruderman (2009). 

The fractional jump in v at the time of the glitch was very large, about 15% (Table 2), 
much larger than the value of 1.7% found by Lyne et al. (2000) to represent most pulsars. 
For PSR B2334+61, Az>c,/z> is dominated by the fastest decay term, with the relative contri- 
butions for the 21-d decay, the 147-d decay and the permanent jump being 85%, 8% and 7%, 
respectively. In most glitch theories, Az>c,/z> represents the moment of inertia fraction which 
is in a superfluid state not tightly coupled to the neutron-star crust (see, e.g., Alpar et al. 
1989, 1993). The large value for PSR B2334+61 suggests that a considerable fraction of 
the star is weakly coupled superfluid, with most of it returning to an equilibrium rotational 
state within about 20 days. About 1% of the rotational inertia appears to be more loosely 
coupled, leading to the 147-d decay, and a further 1% barely coupled at all. 

The very significant quasi-permanent change in v associated with the PSR B2334+61 
glitch (Figure 1) is unusual. Such clear quasi-permanent changes in z> have previously only 
been observed in the Crab pulsar (Lyne et al. 1993). Although it is fairly common for only 
part of Az>g to recover exponentially after a large glitch, in most cases there is a continuing 
approximately linear recovery of v. As discussed above, this was seen in the Vela pulsar 
period, and other examples have been presented and discussed by Lyne et al. (2000) and 
Yuan et al. (2010). For PSR B2334+61, the fractional change Az>p/z> is about 1.1%, about 25 
times larger than the value of ^ 4 x lO"*^ observed for the Crab pulsar. These sudden changes 
in z> implies equally sudden changes in the effective braking torque at the time of the glitch 
since a change in the stellar moment of inertia of this magnitude is ruled out by the fact that 
^Vp/v is several orders of magnitude smaller. The observed change in z> for PSR B2334+61 
seems too large to be accounted for by exterior magnetic field changes as suggested for the 
Crab pulsar (Ruderman et al. 1998; Ruderman 2009). Furthermore, there was no major 
change in the pulse shape at the time of the glitch, suggesting that magnetospheric changes 
are not responsible. Therefore it is most likely that changes interior to the neutron star are 



responsible. The effective braking torque is the difference between magnetospheric torques 
and internal torques believed to be due to vortex creep (e.g., Alpar et al. 1984). This model 
predicts a linear increase in v over long timescales as is seen in the Vela pulsar (Alpar et al. 
1993). For PSR B2334+61 the observed value of i) is two orders of magnitude smaller than 
for Vela, which is in accordance with the predictions of the vortex-creep model (Alpar & 
Baykal 2006). 

The observed long-term step in i) (Table 2) could be interpreted as indicating a third 
exponential decay with a very long timescale. Large positive values of v are commonly 
attributed to the long-term effect of (sometimes unseen) glitches (e.g., Hobbs et al. 2004; 
Alpar & Baykal 2006). For an exponential decay, lS.v^ ^ i^p = ^^ds/Tds- Since we do 
not know Qds = Az^^a/^i^gi we can just set a limit on r^a < 70 years corresponding to 
Q(Q = 1. In principle, the observed post-glitch value of u (Table 1) could be interpreted as 
A'i^'ds and used to solve for Qds and r^a, but it is inconsistent with this, implying a decay 
timescale — Ai^da/A'z/'t^a of less than one year. Evidently the shorter-term decays are not 
exactly exponential or other noise processes contribute to the observed value of 'z>'. 

The oscillations in the post-glitch timing residuals are a remarkable feature of the timing 
behaviour of PSR B2334+61 (Figure 3). These oscillations have a period very close to one 
year, but in Section 3 we showed that an interpretation in terms of position error or proper 
motion is very improbable. Quasi-periodic oscillations in timing residuals have been seen 
in the Crab pulsar (Lyne et al. 1993) and many other pulsars (Hobbs et al. 2010). In none 
of these cases though are the oscillations clearly associated with a glitch although glitches 
have been mentioned as an excitation mechanism (e.g., Zhou et al. 2004; Timokhin 2007). 
The observed quasi-periodic oscillations have often been interpreted as due to Tkachenko 
oscillations in the superfluid vortex array (Ruderman 1970; Popov 2008). Under certain 
reasonable assumptions, these oscillations have a period given by 

Pt - IJlReP^/^yi (5) 

where Rq is the neutron-star radius in units of 10^ cm and P is the pulsar period in seconds. 
For PSR B2334+61, this gives a period of 1.24Rq yr, remarkably close to the observed period 
of 1 yr. This suggests that Tkachenko oscillations are a viable mechanism for the observed 
oscillations. Further investigation is needed to understand the triggering mechanism and 
whether or not the observed oscillation amplitude can be accounted for. 
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Fig. 1. — Glitch of PSR B2334+61: (a) variations of rotational frequency Au relative to the 
pre-glitch solution, (b) an expanded plot of Au where the mean post-glitch value has been 
subtracted from the post-glitch data, and (c) variations of the frequency first derivative z>. 
The vertical dashed line marks the glitch epoch. The plotted points represent fits of u and 
z> to 5 - 10 adjacent TOAs, with Au being the difference between the fitted value of u and 
the value from the pre-glitch solution (extrapolated after the glitch). 
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Fig. 2. — Timing residuals for PSR B2334+61 relative to the pre-glitch solution (Table 1) 
for data up to one year after the glitch. For the upper plot just one exponential decay term 
was fitted, whereas for the lower plot two exponential decay terms were fitted. For further 
details, see text. The vertical dash line indicates the assumed epoch of glitch. 
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Fig. 3. — Timing residuals for the entire seven-year dataspan after fitting for two expo- 
nential recoveries and a permanent changes in z> and i). The pulsar position and frequency 
parameters were held fixed at values obtained from the pre-glitch solution. 
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Table 1. Pre- and post-glitch timing parameters for PSR B2334+61 



Parameter 


Pre-glitch 




Post-glitch 


Right ascension (J2000) 


23*^ 37™ 


09:762(10) 


Declination (J2000) 


+6r51'01'.'53(7) 


Proper motion in RA dirn (mas yr~^) 




-P 




Proper motion in Dec. (mas yr~^) 




-15^^ 




Epoch of position (MJD) 


53100 




Frequency, p (Hz) 


2.01874888482(2) 


2. 


01869363589(2) 


Frequency derivative, z> (lO"-'-^ s~^) 


-780.7074(5) 




-788.332(3) 


Frequency 2nd derivative, i) (10~^^ s~^) 


3.17(6) 




14.46(9) 


Frequency 3rd derivative, z>' (10~^° s""^) 






-0.37(2) 


Frequency epoch (MJD) 


53100.0 




54521.0 


Data span (MJD) 


52486-53609 




53998-55045 


DM (pc cm-3) 


58.410 




58.410 


No. of TO As 


98 




100 


RMS timing residual (ms) 


0.91 




2.07 


Reduced x^ 


1.04 




4.88 



^Prom Hobbs et al. (2004) 
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Table 2. The glitch parameters 



Parameter 


Value 


Glitch epoch (MJD) 


53615(6) 


^vjv (10-6) 


20.5794(12) 


AZ>g/Z> 


0.156(4) 


Az>p (10-15 s-2) 


-8.684(17) 


Az;^ (10-24 s-3) 


8.4(2) 


Az/rfi (/iHz) 


0.19(3) 


Tdl (d) 


21.4(5) 


Az/d2 (/uHz) 


0.119(4) 


T^(d) 


147(2) 


Q = (Ai/rfi + Ayd2)/Az/g 


0.00751(5) 


Data span (MJD) 


52486-55045 


Rms timing residual (ms) 


1.72 


Reduced x^ 


3.32 



